import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 定义区室体积
V_plasma = 5  # 充分灌注区室体积
V_fat = 10    # 脂肪区室体积
V_slow = 5    # 不充分灌注区室体积
V_liver = 2   # 肝脏区室体积

# 定义血流量
Q_plasma = 1  # 充分灌注区室血流量
Q_fat = 0.5    # 脂肪区室血流量
Q_slow = 0.3   # 不充分灌注区室血流量
Q_liver = 0.8  # 肝脏区室血流量

# 定义转移速率
k_plasma_to_fat = 0.1  # 充分灌注区室到脂肪区室的转移速率
k_plasma_to_slow = 0.1 # 充分灌注区室到不充分灌注区室的转移速率
k_plasma_to_liver = 0.2 # 充分灌注区室到肝脏区室的转移速率
k_fat_to_plasma = 0.1   # 脂肪区室到充分灌注区室的转移速率
k_slow_to_plasma = 0.1  # 不充分灌注区室到充分灌注区室的转移速率
k_liver_to_plasma = 0.2 # 肝脏区室到充分灌注区室的转移速率

# 竞争抑制结合参数
Ki_benzene = 1.0  # 苯的抑制常数
Ki_toluene = 0.8  # 甲苯的抑制常数
Ki_xylene = 0.6   # 二甲苯的抑制常数

def pbpk_model(t, y, C_TWA_benzene, C_TWA_toluene, C_TWA_xylene):
    C_plasma_benzene, C_fat_benzene, C_slow_benzene, C_liver_benzene, \
    C_plasma_toluene, C_fat_toluene, C_slow_toluene, C_liver_toluene, \
    C_plasma_xylene, C_fat_xylene, C_slow_xylene, C_liver_xylene = y

    # 苯的浓度变化
    dC_plasma_benzene_dt = (Q_fat * (C_fat_benzene - C_plasma_benzene) / V_plasma +
                            Q_slow * (C_slow_benzene - C_plasma_benzene) / V_plasma +
                            Q_liver * (C_liver_benzene - C_plasma_benzene) / V_plasma +
                            k_fat_to_plasma * C_fat_benzene +
                            k_slow_to_plasma * C_slow_benzene +
                            k_liver_to_plasma * C_liver_benzene -
                            (k_plasma_to_fat + k_plasma_to_slow + k_plasma_to_liver) * C_plasma_benzene +
                            C_TWA_benzene / (1 + C_plasma_toluene / Ki_toluene + C_plasma_xylene / Ki_xylene))

    dC_fat_benzene_dt = k_plasma_to_fat * C_plasma_benzene - k_fat_to_plasma * C_fat_benzene
    dC_slow_benzene_dt = k_plasma_to_slow * C_plasma_benzene - k_slow_to_plasma * C_slow_benzene
    dC_liver_benzene_dt = k_plasma_to_liver * C_plasma_benzene - k_liver_to_plasma * C_liver_benzene

    # 甲苯的浓度变化
    dC_plasma_toluene_dt = (Q_fat * (C_fat_toluene - C_plasma_toluene) / V_plasma +
                            Q_slow * (C_slow_toluene - C_plasma_toluene) / V_plasma +
                            Q_liver * (C_liver_toluene - C_plasma_toluene) / V_plasma +
                            k_fat_to_plasma * C_fat_toluene +
                            k_slow_to_plasma * C_slow_toluene +
                            k_liver_to_plasma * C_liver_toluene -
                            (k_plasma_to_fat + k_plasma_to_slow + k_plasma_to_liver) * C_plasma_toluene +
                            C_TWA_toluene / (1 + C_plasma_benzene / Ki_benzene + C_plasma_xylene / Ki_xylene))

    dC_fat_toluene_dt = k_plasma_to_fat * C_plasma_toluene - k_fat_to_plasma * C_fat_toluene
    dC_slow_toluene_dt = k_plasma_to_slow * C_plasma_toluene - k_slow_to_plasma * C_slow_toluene
    dC_liver_toluene_dt = k_plasma_to_liver * C_plasma_toluene - k_liver_to_plasma * C_liver_toluene

    # 二甲苯的浓度变化
    dC_plasma_xylene_dt = (Q_fat * (C_fat_xylene - C_plasma_xylene) / V_plasma +
                           Q_slow * (C_slow_xylene - C_plasma_xylene) / V_plasma +
                           Q_liver * (C_liver_xylene - C_plasma_xylene) / V_plasma +
                           k_fat_to_plasma * C_fat_xylene +
                           k_slow_to_plasma * C_slow_xylene +
                           k_liver_to_plasma * C_liver_xylene -
                           (k_plasma_to_fat + k_plasma_to_slow + k_plasma_to_liver) * C_plasma_xylene +
                           C_TWA_xylene / (1 + C_plasma_benzene / Ki_benzene + C_plasma_toluene / Ki_toluene))

    dC_fat_xylene_dt = k_plasma_to_fat * C_plasma_xylene - k_fat_to_plasma * C_fat_xylene
    dC_slow_xylene_dt = k_plasma_to_slow * C_plasma_xylene - k_slow_to_plasma * C_slow_xylene
    dC_liver_xylene_dt = k_plasma_to_liver * C_plasma_xylene - k_liver_to_plasma * C_liver_xylene

    return [dC_plasma_benzene_dt, dC_fat_benzene_dt, dC_slow_benzene_dt, dC_liver_benzene_dt,
            dC_plasma_toluene_dt, dC_fat_toluene_dt, dC_slow_toluene_dt, dC_liver_toluene_dt,
            dC_plasma_xylene_dt, dC_fat_xylene_dt, dC_slow_xylene_dt, dC_liver_xylene_dt]


# 初始条件
y0 = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

# 监测数据
C_TWA_benzene = 1.0  # 苯的C_TWA浓度
C_TWA_toluene = 0.8  # 甲苯的C_TWA浓度
C_TWA_xylene = 0.6   # 二甲苯的C_TWA浓度

# 求解ODE
sol = solve_ivp(lambda t, y: pbpk_model(t, y, C_TWA_benzene, C_TWA_toluene, C_TWA_xylene), [0, 8], y0, method='RK45')

# 输出结果
print(sol.t)
print(sol.y)

# 可视化结果
plt.plot(sol.t, sol.y[0], label='plasma benzene')
plt.plot(sol.t, sol.y[1], label='fat benzene')
plt.plot(sol.t, sol.y[2], label='slow benzene')
plt.plot(sol.t, sol.y[3], label='liver benzene')
plt.plot(sol.t, sol.y[4], label='plasma toluene')
plt.plot(sol.t, sol.y[5], label='fat toluene')
plt.plot(sol.t, sol.y[6], label='slow toluene')
plt.plot(sol.t, sol.y[7], label='liver toluene')
plt.xlabel('Time (h)')
plt.ylabel('Concentration (mg/L)')
plt.legend()
plt.show()
